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I. INTRODUCTION 


In recent years functional continuum methods have been extensively used to study the low 
energy sector and the phase diagram of quantum chromodynamics (QCD). These methods 
include functional renormalization group (FRG) flow equations JT]-|3j, Dyson-Schwinger equa¬ 
tions (DSE) |H |5l and variational methods [6]-{8]|- The gauge-variance of the Green's functions 
makes it necessary to fix a gauge, and most techniques such as FRG and DSE initially con¬ 
centrated their studies on the case of covariant gauges. This choice has a two-fold advantage: 
On the one hand, the BRST symmetry and the ensuing Slavnov-Taylor identities provide con¬ 
straints to guide the analysis. More importantly, however, the Kuga-Ojima criterion [9] |T0| 
claims a direct connection, based on the BRST mechanism, between the propagators in Lan¬ 
dau gauge and physical phenomena such as colour confinement. 

At the quantum level, it is not immediately clear if BRST symmetry is naively maintained 
or visible. Most of the functional studies quoted above initially found an infrared vanishing, 
scaling type of solution for the gluon propagator (as the Kugo-Ojima criterion would suggest), 
which is, however, at odds with high-precision lattice simulations fiTUMl . It was later shown 
that infrared finite decoupling solutions could also be obtained, if the infrared behaviour is 
sufficiently constrained fl5l ; such solutions had also been found earlier |[T6V(3 3 . The decou¬ 
pling solutions agree very well with lattice data, but indicate a (soft) BRST breaking in the full 
theory. 

In Ref. HI, we proposed a variational approach that is based on the effective action for the 
gluon propagator. The technique offers several conceptual advantages: For instance, it auto¬ 
matically yields a closed set of integral equations that can be renormalized through conven¬ 
tional counter terms without further truncation. In addition, the variational approach allows, 
in principle, to discern between the scaling and decoupling type of solution, as the one with 
the lower effective action must be realized. Numerically, the approach yields excellent agree¬ 
ment with lattice data that is on par with the best functional methods listed above. 

On the other hand, the variational approach violates BRST symmetry and the Kugo-Ojima 
criterion does not apply so that the question of colour confinment in the variational approach 
is mute at the moment. This explicit violation of BRST symmetry is unavoidable (to a certain 
extent) in any analytical approach, and it is important that the breaking occurs in a controlled 
way. For the variational approach, BRST symmetry is implemented exactly in the target func- 
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tional and an unconstrained variation will give the exact solution. It is only through restric¬ 
tions on the trial variation measure that the breaking of BRST symmetry occurs. As we enlarge 
the variational ansatz space, we improve the quality of the dynamics, and at the same time 
reduce the violation of BRST symmetry (even though we do not have a reliable measure to 
quantify this) Q The situation is, in fact, analogous to the Hamiltonian formulation in Coulomb 
gauge 0 1331436 1 where Gauss' law can be implemented exactly in the target Hamiltonian, but 
only the exact solution of the functional Schrodinger equation will obey Gauss' law exactly. 
Numerically, the variational solution of Ref. HI describes the propagators very well, particu¬ 
larly in the mid-momentum region which is most important phenomenologically, expecially 
for the deconfinement phase transition. It is therefore natural to extend this approach to finite 
temperatures and finite chemichal potentials, and eventually include dynamical fermions. In 
the present paper, we describe the first step in this program, viz. the introduction of non-zero 
temperature. Naturally, this extension has also been done in the FRG l37l and DSE ||38l as 
well as in a perturbative approach Il39l l4QH . The lattice data PTI show a cross-over type of 
signal in the ratio of the two inequivalent Lorentz structures for the gluon propagator, but 
no clear qualitative change at the deconfinement phase transition. In the present paper, we 
will concentrate on the finite-temperature propagators and search for signals of a phase tran¬ 
sition that can be derived from them alone. Eventually, the present work should be extended 
to measure the effective action for the Polyakov loop, which is the real order parameter for 
Yang-Mills theory at finite temperature. This study will be reported elsewhere. 

The paper is organised as follows: In the next section, we revisit the covariant variational prin¬ 
ciple and extend it to finite temperatures through the imaginary time formalism. In section 

the unrenormalized gap equation. Section [Iv] simplifies and improves the renormalization at 
zero temperature discussed in Ref. I§j|. We demonstrate that the same counter terms are also 
sufficient to render the theory finite at any non-zero temperature, and present the fully renor¬ 
malized finite-temperature corrections to the integral equations. In section [VJ we report our 
numerical treatment and present solutions for the ghost and gluon propagator. In particular, 
we study possible signals for a phase transition in the gluon and ghost propagator. In the last 
section, we conclude this study with a brief summary and an outlook on further extensions of 

1 This is consistent since symmetries are invariances of the dynamics, i.e. it is not necessary nor appropriate to 
enforce a symmetry of the exact model onto the truncated dynamics. 


Ill we compute the effective action for the gluon propgagator at finite temperature and derive 
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the variational method. 


II. THE VARIATIONAL PRINCIPLE 


A. The variational method at finite temperatures 


In Ref. 0, we described the basics of the covariant variational principle in quantum field 
theory: for a theory with a quantum field A defined by an action S(A) in Euclidean space 
time, the variation is over all normalised path integral (probability) measures dp(A) used to 
compute quantum averages according to (■ ■ ■ ) f , = f dy(A) ■ ■ •. If the measure is written in 
Radon-Nikodym form du(A) = dA p(A) with a suitable density p, the entropy 

W(/0 = = - f dAp(A) In p(A) (1) 

describes the accessible field space for quantum fluctuations. In the full theory, the entropy of 
the fluctuations is balanced against the classical Euclidean action such that the free action 

F(p) = (S) ft - KW(p) = min (2) 

is minimized. This is because the unique solution of the variational principle (|2]) is the Gibbs- 
like measure 


dyo(A) = Z 1 exp — H 1 S(A) 


dA 


(3) 


whose moments are the conventional Schwinger functions; the minimal free action is simply 
F( m ) — — InZ. The quantum effective action is a constrained version of the free action, 

r(cu) = mmF(y,to) = min|(S) F - KW(}i) \ (Q) fi =o;|, (4) 

where an arbitrary operator Q is fixed to a prescribed external value to. Usually, Q is chosen 
as the quantum field operator A itself whence to becomes the classical field, but this is not 
mandatory: for the present study, it is more convenient to take Q as the 2-point function of 
the quantum field and compute the optimal propagator from the overall minimisation of the 
effective action, 5Y/5to = 0. 

To transcribe this method to non-zero temperatures, we can follow the standard imaginary 
time formalism and restrict the Euclidean time interval to [0, f>] with anti/periodic temporal 
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boundary conditions for fermions/bosons, respectively. This procedure works because the 
free action in our approach is always computed from a conventional quantum field theory, in 
which temperature can be introduced as described. The minimal free action eq. |2]) obtained 
for the Gibbs measure is still F(/to) = — InZ, but with the partition function 


Z(j8) = 


dA exp 


—fT 1 S(A) 


(5) 


now computed using fields with the appropriate temporal boundary conditions. With the 
corresponding Fourier decomposition (in the bosonic case) 

A(t,x) = r 1 E /7S3 ei(v ” t+kx) MV, Vn = 2 a n (nez) ( 6 ) 

nez j \ Ln ) P 

and translational invariance (in both space and time separately), it is easy to see that the space- 
time volume factorizes, = — In Z(j8) = fiV • /. The free action density / obtained in this 

way agrees, in the free case, with the thermodynamical free energy density of a non-interacting 
grand-canonical Bose gas (including the self-energy), provided that the path integral measure 
was properly normalize^ H42l . This observation generalizes to interacting field theories, and 
in particular to gauge theories. In the latter case, the inclusion of the Faddeev-Popov (FP) 
determinant is crucial, and the Matsubara frequencies for the FP ghost fields must involve 
even multiples of n / f J > (as in eq. Q) even though the ghost fields obey Fermi statistics C42l . 
The inclusion of the ghost fields also requires an adaption of the variational principle: since 
the FP determinant J (A) describes the natural measure on the orbit space of the gauge fixed 
theory, the entropy of all trial measures must be computed relative to the FP determinant, 
which amounts to replacing the entropy by the relative entropy 0 


W{}i) = W^i) + (ln(.7)) fi = - (In ip/J)) H . (7) 


B. Gaussian trial measures 


The zero-temperature investigations in Ref. 0 have shown that a Gaussian ansatz for the path 
integral measure d}i is sufficient to describe the propagators of the theory accurately. We will 
therefore try a similar ansatz at finite temperatures. 


dji(A) = A f -dAJiA) 


1—2ft' 


• exp 


,(*/ y) K (*) 5 ab & v (*, y ) A v iy) 


( 8 ) 


2 This equality holds in units where ft = 1 and the Boltzmann constant k B = 1, which we assume from here on. 
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where a is a variational parameter and the normalisation Af depends, in general, on a and on 
the kernel to. We have also assumed that the kernel to can be taken colour diagonal because 
of global colour invariance left after (covariant) gauge fixing. Finally, we have also introduced 
the abbreviation 



(9) 


for the spacetime integral at finite temperature. The corresponding momentum integral is 




( 10 ) 


i.e. the integral over /cq is always understood as a discrete sum ove the Matsubara frequencies 
ko = v n = 2nn/f>. The overall physical picture conveyed by the ansatz eq. (|8]> is a weakly 
interacting (constituent) gluon with an enhanced weight (for tx. > 0) near the Gribov horizon. 


It is clear that the gauge field must be Fourier decomposed as in eq. j6); similarly, the FP 
determinant for covariant gauges 


J(A) = Det 


ab 


-d u D 


/Det 


-DS‘ 


ab 


= Det 


-ns ab -gd F f 


abc 


A \ * 


/Det 


-□<? 


ab 


( 11 ) 


must be evaluated with periodic boundary conditions in time (even though the ghosts obey 
Fermi statistics). 


Let us next look at the gluon propagator 

< A;(x) A b v (y) ) = [ dke ik (*-*) • 6 ab D }lv (k) . (12) 

At finite temperatures, D fiV (k) is still a symmetric rank-2 tensor, but the overall 0(4) Lorentz 
invariance is broken because the heat bath singles out a rest frame specified by u }l = (1,0,0,0). 
Assuming that the remaining spatial 0(3) symmetry remains unbroken, D }iV (k ) must be a lin¬ 
ear combination of all 0(3)-invariant symmetric rank-2 tensors that can be formed from k fl 
and u fl . Furthermore, BRST invariance entails the Ward identity for the longitudinal propaga¬ 
tor. 


kW l D flv (k) =£, 


(13) 
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where £ is the coefficient for the gauge fixing term in covariant gauges, £g x = X (9 ;/ A") 2 . 
From all these constraints, the gluon propagator has just two remaining Lorentz structures]^] 

D, w (k) = D ± (k ) P±(k) + D„ (k) p| v (*) + Y2 kj w (14) 

D^(fc) = DJ 1 ^) 7^(*) + D^(fc) ?]!„(*) + ^ ■ (15) 

The two projectors V^ and V are both 4-dimensionally transversal, but V is also 3- 
dimensionally transversal, while V is 3-dimensionally longitudinal, 

V^ik) = (1 - S }l0 ) (1 - S v0 ) (V - (16) 

V\Uk) = Uv - k -p-) - V^ v (k 0 ,k ). (17) 

These projectors enjoy the usual properties (V L ) 2 = V L and (pH ) 2 = pH, as well as pH • V 1 = 
P- ■ V =0; the traces in four euclidean dimensions are tr V L = 2 and tr V = 1 . Obviously, 
V L + pH = V is the usual 4-transversal projector. 

Returning to the Gaussian trial measure eq. | 8 |, we note that, for a = the kernel to equals 
the inverse gluon propagator and is thus subject to a Lorentz decomposition similar to eq. ( [T5| . 
For ol X \r the FP determinant J contributes; below, we will treat J in the so-called curvature 
approximation Il35l which maintains the Lorentz structure of the kernel. As a consequence, 
the kernel must have the general form for all a, 

<x> }lv (k) = co±(k) -V^yik) + cu\\(k) -Pjv(fc) +j~ 1 k F k v , (18) 

where 7 = £ up to radiative corrections. The two scalar dressing functions to_ and to depend 
on kjt only through the two invariants ko and |k|. 

C. Curvature approximation and the ghost DSE 

To complete our computational tools, we must also evaluate the FP determinant J{A) from 
eq. ( jTT] l and its expectation value in the trial measure (| 8 j. Unfortunately, this cannot be com¬ 
puted in closed form and we will resort to the finite-temperature equivalent of the so-called 

3 To simplify the notation, we will often write f(k) instead of /(fco,k) for a general function of the 4-momentum 
k a = (fco,k). If a function is 0(4) symmetric, f(k) is really a scalar function of k 2 = Ajj + k 2 ; the distinction 
should always be clear from the context. 
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curvature approximation which was previously shown to be effective both in the Hamiltonian 
approach to Coulomb gauge f35il . and in covariant gauges at zero temepratures f6j. In this 
approximation, the FP determinant is expressed as 


In J[A\ ~ f d{x,y) x%(x,y) ■ A^x) A%), (19) 

where the curvature is given by the expectation value 

^ ) (20) 
taken with the Gaussian measure eq. (|8jl. From global colour invariance, we have again — 
5 nb X}iv- Furthermore, the arguments put forward above to determine the Lorentz structure of 
the propagator also apply to the curvature. As a consequence, we have 

xMk) = X ±(k) ■ V^k)+ Xll (k) • Vf n ,(k) + • • • , ( 21 ) 

where the dots indicate irrelevant 4-longitudinal pieces which vanish in Landau gauge. We 
will derive integral equations for the profile functions in the next section. 


If we put everything together, our Ansatz |[8jl for the trial measure becomes 

~\J dkA;(-k) {co ± (k) • V^,(k) + d)|| (k) • v\ v (k) + T (k) • k }l k v } 
where the profile functions are now 


dp(A) = Af ■ exp 


, ( 22 ) 


<Wl = + (1 - 2a) xl 

®ll = w ll + (l - 2 a )*I* ^ 


A few remarks are in order at this point. Firstly, the full gluon propagator in our approach is 


D$(k) = S' 


ab 


(k) 1 • VL (k) + d>ii (k) 1 • v\\ v (k) + r (k) 


k }l k v 


}IV 


k 4 


(24) 


The 4-longitudinal (last) piece in this expression is not very important and will even vanish in 
Landau gauge £—>■(), which we will study exclusively below. 

Secondly, the operator held fixed in our trial measure is not the gauge field itself, but rather 
the gluon propagator (cf. eq. p4|), i.e. according to the general setup explained in section 
our variation principle optimizes the effective action for the gluon propagator. This is appropriate 
as long as we are not interested in higher vertices. 


IIA 
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Finally, one might think that the curvature is irrelevant, as the final form of the trial measure 
including the FP determinant is still Gaussian. This is not so because the curvature reappears 
in the relative entropy W that enters the free action. 


D. The ghost sector 


The ghost propagator is the expectation value of the inverse Faddeev-Popov operator which 
we decompose according to G 1 = — D /( = G (] 1 + h, where 

[Gq- 1 ] ab = -□ 5 ab , h ah = gf aic d fl A c }l (25) 

with the structure constants f abc of the colour group SU(N). From this, it is easy to see 
that the ghost propgator (G) satisfies the Dyson equation (G)“ 1 = G 0 1 + (hG) (G) _1 . The 
expectation value (hG) involves the ghost-gluon Green's function, which in turn is obtained 
by attaching full propagators to the proper ghost-gluon vertex l abc . Using roman digits for 
the combination of Lorentz, colour and space-time indices, the ghost gluon vertex is defined 
by 

( A( 1 ) G(2,3) > = —D(l, 1 ') < G( 2 , 2 ') ) r(2',3' ; 1 ') ( G(3'3) >. (26) 


If we also introduce the ghost form factor ij(k) as the deviation of the full ghost propagator 
from the free one. 


(G) ah (x,y) = I dke ik 5 ab ^, 


(27) 


(28) 


we find the exact momentum space relation 

s" b n (k)- 1 = S‘ b +gf bi f d, ^ D Zf-pA r"(, -k,k). 

Here, we have used global colour invariance to deduce the colour structure = <5 ab D flv for 
the full gluon propagator, and the momentum routing in the vertex is 


F t cb (c,-k,k) : 



A a 


9 




The exact eq. ( [28) is referred to as the ghost Dyson-Schwinger equation (DSE). In the rainbow- 
ladder approximation, the proper ghost-gluon vertex is replaced by the bare vertex defined by 


r 0 (2,3;l) 


SG- 1 (2,3) 

M(l) 


(29) 


which reads, with the momentum routing above. 


[ r o ] d ' b (q ~k,k)^ igf dcb (q - k) v . 


(30) 


The ladder approximation has been proven to be very reliable at T = 0 and we expect it 


(31) 


to hold equally well at non-zero temperatures. From eq. (281, we then get a closed integral 
equation for the ghost form factor, 

*(*)-' = 1 - Ng 2 l d-l 

To complete the derivation of the ghost sector, we have to insert the gluon propagator eq. 
from our Gaussian ansatz into eq. ( [31) . This simplifies considerably in pure Landau gauge 
(£ = 0) when the gauge connection and the gluon propagator are 4-dimensionally transversal. 
From this point on, zve will therefore study Landau gauge (f = 0) exclusively. In this case, the 
ghost DSE ( |3T) becomes 


p ( k-q) 2 <D\\(q) 

at 2 k2 f * h( k ~ d) 
-Ng>- ¥ t q ' 


(l-(k-q) 2 ) hx(i|) 1 — a)|| (q) 


-1 


(32) 


(k-ci ) 2 

where k 2 = Icq + k 2 is the 4-momentum square and we have also introduced four- and three 
dimensional unit vectors 


v — _ 

}l ~W 2 ' 


k- — 

Ikl- 


(33) 


Next, we concentrate on the curvature eq. 


* (1 - 2 > = Tr '< G ^- G 


. In our compact roman digits notation, we have 
S(-dD) 


M(l) 


The derivatives yield the bare ghost-gluon vertex To ( [29) , which is field-independent. To the 
given (formal) loop order, we therefore have (C To G) ~ (G) To (G). The remaining expecta¬ 
tion values are merely the ghost propagators computed earlier. Restoring all arguments, we 
obtain the momentum space expression 


xU k ) = -* ab 



f l(k-q )q(q) 
(k-ciYci 2 


(k ~ q)}i qv 


(34) 
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In Landau gauge, only the 4-transversal parts of this tensor contribute, which we find by 
contracting with the corresponding projectors V and VK respectively After some straight¬ 
forward algebra, the profile functions in eq. pT) are 


y , (k ) _ 1 N? 2 [ da ?(*-</) */(</) q 2 

X±W- 2 N§ jdq (k _ q)2 ^ 


X\\ ( k ) = N g 2 I d£ J 


( k~q ) 2 
(k- q r 


!- (k-q)^ 


1 


1- (k-q) 


(35) 


III. THE GAP EQUATION AT FINITE TEMPERATURES 

A. The classical action 

To derive the effective action for the gluon prop gator, we must first evaluate the expectation 
value of the gauge-fixed YM action in the trial measure eq. (|8j. 



The relevant correlators are easily evaluated using Wick's theorem, global colour invariance, 
and the SU(N) relation f ahc f nbc = N(N 2 — 1) to perform the colour traces. In momentum 
space, the result takes the form 


( S g f ) (N 2 — 1) [5V j dk (k 2 5 }W - (1 - r 1 ) k n k v ) D fW (k) 


(37) 


+ ^(N 2 -1 )flV 


dk D m {k) 


Ng 2 


(N 2 — 1) (5V d(k,q)D nv (k)D V}l {q), 
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where V is the 3-dimensional space volume. Next, we insert the representation (241 and take 
the Landau gauge limit £ —> 0 to fincj^ 


( Sgf) 


= - / dkk 2 


(N 2 - 1) p>V 2 Jp 


d-1 1 

+ 


co±(k) oo\\ (k) 


Ng 2 


dk 


d-1 1 

+ 


co±(k) co\\(k) 


¥ / d <M){ 1+ ( , £ : y, + _ 4 ,, 

4 Jfl u \co±(k)u>±(q) co±(k)co\\(q) 


d — 2 — (k • q) 


1_3_ 

<? 2 


+ 


cum (k) to || (q) 


(fc-if) 2 + (l - (k • q) : 


k 2 q 2 

-1 4-1- -L- 

+ k 2 + q 2 


(38) 


To further simplify this expression, we can exploit the remaining rotational symmetry to 
rewrite some of the integrals in which the integrand depends on 0(3) invariants only Af¬ 
ter some straightforward but lengthy algebra, the classical action can be recast to 


(S gf ) =^(N 2 -1)/3V Jdqq 2 


0 ±(<?) On (q) 


+ 7 Ng 2 (N 2 - 1) §iV [ d(k,q)f _ A . + _ + _ ^'_ q \ }, 

4 Jp \co±{k) ui ± {q) co±(k) cuy (q) to\\ (k) u)\\ (q) j 


(39) 


where the coefficients are given, for d = 3 space dimensions, by 

-I 




C(k, q) — - + 


2 1 /k 2 ~ 2 


3 3 U 2 


+ 


4 k 2 q 2 
3 fc 2 q 2 


(40) 


B. The entropy 


As the last ingredient, we need the relative entropy of the path integral measure ( [22} with 
respect to the Faddeev-Popov determinant. From eq. (j7| and the curvature approximation 


4 We have dropped the field- and temperature-independent constant 


1 n r V 

2 fiV (N 2 — 1 ) Jdk=- 


(N 2 -l) E / 

neZ J 


d 3 k 

(27l) 3 ' 
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eq. (JT9JI, we obtain 


W = (-]np) + (]nj) 

« -lnA/'A^ j^d(x,y)d^(x,y)(A a tl (x)A b v (y))-^ J^d(x,y) xf v (x,y) (A a F (x) A b v (y)) 
= - InA/" + \J d(x,y) (A a F (x)A b (y)) {d>J (x,y) - xft v (x,y)} 

= - Indet + * (N 2 - 1 ) J d(x,y) G)^(x,y) {dy (x,y) - ;fyv(*,y)} 

=_ 2 Trin ( 27 t) + 2^ N2_i ^ y 


where the explicit form A/” = det[d>/ (2zr)] 2 for the normalisation in eq. (22 1 was used. In pure 
Landau gauge in three space dimensions, we have o) ;a ] (k ) dj fiV (k) = tr P(k) = 3, where V is 
the 4-transversal projector. Likewise, 


= ti A 1 {k)x±(k)-trV ± (k) + w ll 1 (k) x\\ (k) • tr 7? 11 (fc) 

= 2d> ± (k)- 1 x±(k) +d>||(k)' 1 2C||(k). 


After dropping an irrelevant temperature- and field-independent term, the entropy becomes 


W = 



1) Tr In 


2n > lv 


^l p \\ 

2n ' ' 


}IV 


i(N 2 -l) pvjM{2 


X±{ k ) 

&±(k) 


*11 W f 

cum (k) J ’ 


Since the projectors P~ and 'P are orthogonal, the determinant in the first term factorizes. 
The dimensions of the subspaces on which P- and V project are 2 and 1, respectively, and 
in each subspace, the projector acts as unity. Expressing the functional trace in momentum 
space, we finally obtain 


W = --(N 2 - l)pV ( ct/ci 2 ln^^ +ln^-^ 

2 k r Js I 2tc 2 tc 


x{ k ) , d ( k ) 


co±(k) co\\(k) 


(41) 


The difference of eqs. (39' and (41 1 is the effective action T((v±,co\\) for the (inverse) gluon 


propagator within the curvature approximation ([19). 
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C. The gap equation 


It is now straightforward to compute the optimal kernels (or gluon propagator) from the gap 
equations 


<5r(cuj_,a;||) _ 6T(u)±,a)\\) _ 
5co ±(k) to||(k) 

After more algebra, we obtain for the variation w.r.t co_ (k) 


(42) 


0 = k 2 + ^Ng 2 


A 


^ B(q) 
p q [&±(q) + &±{q)_ 


- (V ± (k) + x±(k) + f dq 
H 




p 1 cvj_(q) Scv^k)- 1 ' 


(43) 


where the coefficients A and B(q) are given in eq. (40'. The curvature depends on co±(q) 
through the ghost form factor rj which enters eq. ( J35) , but this is a higher loop effect that can 
safely be neglected within the curvature approximation. Dropping the last term in eq. 
we arrive at 


u> ± (fc) = kl + k 2 + x± (k) + M\ (p) 

A , B(q) 


M 2 (/3) = ^ Ng 2 j dq 


d’Aq) ®|| (q) 

The second gap equation is derived in exactly the same fashion and takes the form 

B(k) C(k,q) 


(44) 

(45) 


d>|| (k) = k 2 0 + k 2 + x\\ (k) + Ng 2 J dq 


®±(?) 6>\\(q) 


(46) 


In constrast to the previous case (45I, the last piece has no direct interpretation as a mass term 
as it depends on the external momentum k ?t . However, at k = 0 (with k o 7 ^ 0), we find 


Ng 2 / dq 


B(k) , C(k,q) 


k =0 


Mi 03), 


Ip [ a>±(q) cv ± (q)_ 

as can be shown by explicit calculation. We can thus factorize the momentum dependence in 
the last term on the rhs of eq. |46) and recast the second gap equation to 


k 2 


d>\\ ( k ) =k 0 + k-+ X\\ (k) + Mj_ (/3) + k2 - k2 AM (| (P) 


A Mp) = 1 -Ng 2 j dq 


d>±(q) 


0 

^~3q 2 
flo + q 2 


a w (q) 


(47) 


(48) 


, (441 and (W7) form a closed set of integral equations to determine the gluon and 


Eqs. 

ghost propagators in Landau gauge at finite temperatures. 
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IV. RENORMALIZATION 


Next we turn to the renormalization of our integral equation system. We expect that (i) 
the equations reduce to the known zero-temperature case in the limit (5 —> oo and (ii) the 
counterterms fixed in the zero-temperature limit are sufficient to renormalize the system at 
any temperature. 

A. The zero temperature limit 

Let us first study how the zero temperature limit formally arises in our integral equation 
system Q The sum over Matsubara frequencies can generally be rewritten by Poisson resum¬ 


mation. 



(49) 


By the Riemann-Lebesgue lemma, the Fourier transform fix) vanishes at large arguments 
|x| —» oo, i.e. only the term with m = 0 contributes in the zero-temperature limit. 



(50) 


Let us assume for the moment that the ghost form factor //(7c) approaches the 0(4)-invariant 


zero-temperature form factor t]o(k), cf. appendix A1 (This assumption will be justified a 


posteriori below). From eq. (j34j) and the limit (50), we have 



Contracting with V^ v (k) gives 



where Xo is the zero-temperature curvature, cf. appendix |A 1| From the second equation (35), 


it follows that x\\ (k) y 3 Xo(k) — 2 lim x±(k) = Xo(k) in the same limit and we have 

P. _ 


lim x±(k) = lim x\\(k) = Xo(k) 

fi —>oo fi—>oo 11 


(51) 


5 Since we are dealing with the unrenormalized equations, an 0(3) invariant cutoff is implicitly assumed in this 
section. 
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In appendix A 2 it is further shown that the mass functions in the gap equation have the 


low-temperature limit M]_(/3) —* Mg and AM»(|3) —> 0. From eqs. (441, (47 1 and (51i, it then 
follows immediately that 


lim to ± (k) = lim to n (k) = kg + k 2 + Xo (k) + Mg = to 0 (k) 

ft— >oo /?—>00 11 


(52) 


Finally, this relation shows that the second line in eq. (321 vanishes as ft —? oo, while the first 


line approaches ?/o (k) as can be seen from eqs. (A1' and (A2 ■. Thus, the self-consistent solution 
has also 


lim tj(k) = 7/o (fc), 

/3->oo 


(53) 


which justifies our initial assumption above. Alltogether, eqs. (511-( 53) are the expected ft —? co 
limit in which the full 0(4) symmetry is restored. 


B. Renormalization at zero temperature 


In this section, we recall (and simplify) the renormalization of the zero temperature system 
as layed out in Ref. This information is required later on, since we must fix all counter 
terms at zero temperature in order to compare different temperature settings reliably. After 
introducing ghost fields c,c to make the Faddeev-Popov determinant local, three counter 
terms are required, 

Oct = i (9 V A% - 9 v a ;) 2 + 5M 2 • i (a ;) 2 + SZ C • 9 F cd^c (54) 

corresponding to a gluon wave function, gluon mass, and ghost field renormalisation. (No 
vertex renormalisation is induced by the theory.) The counterterms f d 4 x £ ct must be added 
to the exponent of the (Gaussian) trial measure eq. (|8]> used to compute the 77-point functions, 
i.e. they are not counterterms to the Yang-Mills action directly. The renormalized integral 
equation system at zero temperature then becomes 

c!>o(k) = k 2 + Ml +Xo{k) Mq = Ng 2 iff + 5M\ 

7/0 (fc)- 1 = 1 -Ng 2 if (k) - 3Z C X o(k) = Ng 2 if (k) +S X + k 2 SZ A , (55) 


where the explicit form of the various loop integrals is listed in appendix A1 The mass 
counterterm has been split, SM 2 = 5M\ + d X , between the mass and the curvature since only 


the sum of these contributions enters the gap equation (551. 
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As layed out in Ref. [61, we have to choose slightly unusual renormalization conditions to 
overcome numerical problems in the deep infrared. Normally, one would prescribe the value 
of the propagators and vertices at a common (large) scale g 1, except for the mass countert¬ 
erm which determines the propagator at zero momentum (in Euclidean space). In the present 
case, we impose the value of the ghost propagator, or ghost form factor g{gc), at a low scale 
g c < 1 (which could even be taken to be g c = 0) in order to discern the scaling from the 
decoupling solution [6|; this condition determines the counter term 

5Z C = 1- Ng 2 tf\g c ) - rjoipc)- 1 ■ (56) 


The remaining renormalization conditions for the gluon can be taken conventionally: we im¬ 
pose the value of the (inverse) propagator at a large scale /(> 1 to fix the field renormalisation 
5Za, and at a small scale go <C 1 to fix the mass counterterm^j 

u> 0 (g)=:Zg 2 , d>o(f/ 0 )=:ZM 2 a . 

For later reference, we note the explicit form of the counterterm coefficient^] 


(57) 


u 2 — m 2 a ? 4°4/0 _ 

1 + SZ A = Z t — - £~Ng 2 - x VM x Kt0> 


r(°)/ 




Ho 


H z 


Ho 


Ml + 5M 2 = Z g- 


M a ~ Ho Ar „ 2 H 2 4 0) (Ho) ~ Ho 4 0) (H) 


2 ~ N S 

r - Ho /' z 

which complement eq. |56] l in the ghost sector. With these conditions, the renormalized inte¬ 
gral equation system at zero temperature become^] 


Ho 


(58) 


Vo( k ) 1 = Vo (Ho) 1 ~ Ng 2 4 ° } (k) ~ 4 ° \Hc )] 


(59) 


co 0 (k) = Z 


H 2 ~ m a 


H z 


hI 


A k 2 


z 4 g 


N g 2 




H 2 o 


H 

r(0) / 


Ho 


H 2 (lx\k) W)- 


( 0 ), 


- k 2 (4 U V) - 4 0) (Fo)) - Ho (4V - 4 u> (h)) 


( 0 ), 


r(0). 


Note that the difference of integrals in the square brackets is finite: For the first equation 
this is clear because 1^ is logarithmically UV divergent by power counting. In the second 

6 For the decoupling solution, we could choose go -> 0 so that M 2 A could be interpreted as a constituent gluon 
mass; in the general case, however, we must take go > 0 to avoid infrared singularities, and M 2 A becomes a 
general mass parameter. 

7 The symbol Z should not be confused with the gluon field renormalization constant, which is Z A = 1 + SZ A - 

8 The finite parts in the ghost counter term SZ C are such that the first equation in 59 could also be obtained by 

simply subtracting the bare ghost DSE for t]o(k) at the renormalization scale g c - 
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equation (59 > l' x 0) is quadratically divergent, but the subleading logarithmic divergence is also 


eliminated due to the clever combination of fj°* at three different scales. 


The system eq. (591 is well suited for a numerical evaluation. It contains, besides the three 
renormalization constants {tjo(nc), Z, M 2 a }, also the coupling strengthen the combination 
Ng 2 . However, this parameter is redundant. To see this, consider a solution dy(/c), ijq (k ) of 


eq. (59' obtained with a set of parameters {ijo(}i c ),Z, M\, Ng 2 }. From the form of the loop 
integrals in eq. ( |59| listed in appendix |A 3} it is apparent that the rescaled parameters 

{ f/oOc), Z, M 2 a , Ng 2 } ->■ { a ■ t]o(ii c ), b • Z, M 2 A , b/a 2 • Ng 2 } (60) 

lead to the solution ar]o(k) and bcoo(k), in which the propagators are scaled by arbitrary 
factors. Thus, any change in Ng 2 can be compensated by a corresponding rescaling of the 
propagators (and vice versa), including their field renormalisation constants f]o{pt c ) and Z. We 
can thus choose to set Ng 2 = 1 in eq. (J59]l, with the provision that the propagators (and their 
field renormalization constants) may later be multiplied by an overall constant when fitting to 
lattice data. The solution of eq. ( [59] > therefore depends on the three renormalization constants 
{ tjoigc), Z, M 2 a } only, where the first two determine the overall scale of the propagators, while 
M a determines the behaviour in the intermediate momentum region. Formally, M A can be 
chosen arbitrarily, but it is related through a gap equation to the gluon propagator and the 
bare Mass Mo- Since this gap equation involes Ng 2 , different colour groups SU(N) will 
require, in general, different value for M 2 A . For the case of SU (2), we will give a numerical 


estimate for M A from high-precision lattice data in section V A below. 


C. Renormalization at finite temperatures 


Next we turn to the integral equation system at finite temperatures. In the previous section, 
the counter terms were constructed at zero temperature and it is important to use the same 
counter terms also at non-zero temperatures in order to compare different temperatures reli¬ 


ably. To do so, we merely have to add {—5Z C ) from eq. (561 to the right hand side of eq. 


and likewise ( 5Z A k 2 + SM 2 ) from eq. (58 ' to the right hand side of eqs. (44' and (47'. The re¬ 
sulting system of equations can be written in two equivalent forms. This is best demonstrated 


In the present approximation scheme, Ng 2 is a finite adjustable parameter since there are no vertex corrections. 
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at the ghost equation: 


k 2 


V(k) =1 -Ng Ijj( k) — Ng Tj 


V(k-q) 


k 2 Jfi {k-q) 2 


!- (k-q) 2 ] [d f ± (q) 1 - G\\(q) 1 


-<iZ r 


mLr,(k) 


k 2 


= 1 - Ng 2 I n (k) - Ng 2 ^ - 1 - Ng 2 4 0) (p c ) - J/o(^c) 


fc 2 


k 2 


= Vo(Fc) ~ Ng yl v (k) - P v \y c )\ ~ Ng ^ L,,{k ) 
= >1o{k)- 1 - Ng 2 \l n (k) - 4° } W] - Ng 2 ^ L v (k), 


(61) 


where the zero temperature equation (591 was used in the last step. The advantage of this 
formulation is that it avoids all explicit renormalization factors, as the finite temperature solu¬ 
tion is expressed as its zero temperature counterpart plus a temperature-dependent correction 
which must be finite as all counterterms have already been consumed in j/q (k). To confirm 
this, we note that the curvature and mass term in the gap equations |44| , ( [4 7} are subdomi¬ 
nant, which implies that \co±(k)^ 1 — d>i|(A:)~ 1 | ~ 0(fc~ 4 ) at large momenta and (k) is finite 
by power counting. As for the difference — ljj°'(k) of logarithmically divergent loop inte¬ 
grals, we note that this would be finite at zero temperature. At nonzero temeprature, however, 
this statement is less obvious, since I,](k) and ij^ (7c) involve different form factors in their 
integrand. In principle, the restauration of 0(4) invariance at large momenta should ensure 
that the difference is finite at any temperature, but this limit is very hard to reach numerically. 
We will therefore follow a different route in our numerical investigation below. 

Similar considerations can also be applied to the gap equation. After some algebra, we obtain 


>1 (kor | k |) 1 = qo(k) 1 -Ng 2 |l,(fc 0 ,|k|) -tf\k) 
Q ± (k 0 ,\k\) = G 0 (k) + Ng 2 \l^(k 0 ,\k\) - lf\k) 


k 2 


* 2 o 


bk 2 

M 2 G6)-M 2 


Ng 2 Lfi(ko, | k |) 


(62) 


d>ll(*o,|k|) = a>o(*) + Ng 2 I"(k 0 ,\k\)-lj°\k) + M 2 QS ) - M 2 0 


k 2 


fc 2 + k 2 


A M 2 (/l). 


To avoid confusion, we have reinserted to full momentum dependence in the finite-tempera¬ 
ture profiles and loop integrals. Explicit expressions for the loop integrals entering this system 
can be found in appendix |A3| 

As mentioned above, the system eq. ( |62| is hard to treat numerically because the difference 
of loop integrals involve non-zero and zero temperature form factors, respectively. We will 
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therefore rewrite these equations to bring them as close as possible to the zero-temperature 
case. Let us again consider the ghost equation as an example. From the second line in eq. (|6TJ, 
we have 


y(k 0 , |k|) 1 =r}{Q,n c ) 1 — Ng 2 I t1 (k 0 , |k|) - l v (0,}i c ) -Ng 2 


k 2 


L * 2 o 


k 2 


Li](kor |k) 


Here, we have chosen the renormalization point at vanishing Matsubara frequency^ The 
benifit of this reformulation is that it only involves the difference of the same finite-temperature 
loop integrals (with the same form factors) at different scales. The cumbersome difference of 
loop integrals with different profiles is now hidden in a single number, namely the (finite) 
value of the form factor at the renormalization scale. 


V( 0,Hc) 1= 7/o(jM C ) 1 


i v (o,}i c ) - 4 0) (fc) - Mo ,}ic) 


(63) 


Again, we cannot compute this quantity reliably due to numerical issues, and we cannot 
choose it arbitrarily either, as it is determined completely by the zero-temperature counter 


terms. We will instead fix ij(0,pi c ) indirectly as explained in detail in section VB below. 


The gap equations can be rewritten in much the same way. After some lengthy algebra, we 


10 The orientation of the renormalization point in 3-space is arbitrary due to the residual rotation symmetry, and 
we can choose the same direction as the external momentum, (0, }i c ■ k), if this is ever necessary. A similar 
statement applies to the other two renormalization scales }Iq and }i. 
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arrive at the following system: 


V(ko, k) 1 = // (0 ,}i c ) 1 - Ng 2 I v (k 0 , k ) - I v (0, y c ) - Ng 2 -_ 2 L ?/ (fc 0 , fc) - L v (0, y c ) 


k 2 


cv ± (k 0 ,k) = cv±(0, y) + 


F 2 ~ Fo 


cv ± (0,y 0 ) - co±(0,y) 


+ 


Ng 2 


2 2 t F 

f 2 - n 1 


I (k 0 , k) - I (0, y 0 ) ~(k 0 + k) I x (0, y) - I (0, y 0 ) 


Fo 


#(*„,*)-1* (0,p) 


v\\(ko,k) = O7|i(0.,/<) + 


F 2 -(k 2 0 +k 2 ) r 


f 2 - fI 


kl 


a \\(0'Fo) - ^||(0,/OJ - k 2 + k 2 AM ||(^) + 


+ 


i „ 2 

M 2 -/'o 1 


4(fco,fc)-1,1(0,^) ~(k 2 0 + k 2 ) I"(0,y)-I"(0,yo) 


■II, 


■II, 


Fo 


4m - 4(o, f) 


(64) 


Here, k = |k| and we have spelled out all arguments to the loop integrals for clarity. The gap 
equations now exhibit the same combination of curvature integrals at three different scales 


which eliminated all divergences in the zero-temperature case. Alltogether, eq. (641 contains, 
in addition to the coupling constant Ng 2 , five temperature-dependent constants. 


co±(0,y 0 ), CJ ± (0,}1 0 ), 6l7|| (0 ,}1 0 ), to\\ (0, yo), ij(0,Fc), (65) 


which are all determined, in principle, by the zero-temperature counter terms through finite 
equations that are, however, hard to treat numerically. For completeness, these relations are 
listed in appendix |A 4} although we will not use them and determine the constants in eq. |65} 
indirectly by other means. 


The systems (621 and (|64]> are the main result of this paper. In the next chapter, we will give 


details on its numerical implementation and discuss the solutions. 
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V. NUMERICAL TREATMENT AND RESULTS 


A. Renormalization at zero temperature 


As explained in the last section, the renormalization of our approach must be carried out 


entirely at zero temperature. In section IV B we have presented a simplified procedure, where 


all counter terms are determined by simple local prescriptions to low order Green's funcitons. 
Since this is different from the procedure initially layed out in Ref. |6j, we cannot directly 
take over those results and must first re-determine the optimal renormalization parameters to 
match the known lattice results. This is carried out at zero temperature, and we take again 
the high-precision SU(2) lattice data from Ref. Ifl3l . 

In figure [TJ we compare the best fit of our zero-temperature formulation to the lattice data 
for both the gluon propagator and the ghost form factor. Since we already know that the 
physically realized solution is subcritical (decoupling solution), only this type of solution is 
investigated. The scale }Iq for the gluon mass renormalization could thus be taken to zero, but 
it is numerically favourable to impose this condition at a small but non-zero scale }Iq ~ /i/10. 
As can be seen, the new renormalization prescription describes the lattice data with similar 
(high) precision as the more complicated prescription in Ref. f6|. We find Ma / }t ~ 0.49 at 
}i ~ 1 GeV to be a favourable parameter, but this determination has rather great uncertainties. 
As can be seen from figure [IJ this is not because the approach could not describe the data, but 
rather since different combinations of the renormalization constants produce similar results, 
i.e. there are (almost) flat directions in renormalization space. From the comparision to the 
lattice data, we can also determine our parameters in absolute numbers. For the choice in 
figure [lJ we can match the data at the point }t = 0.96 GeV, where we find 


M a « 0.48GeV. 


( 66 ) 


The renormalization scale }i is, of course, arbitrary as long as it is sufficently deep in the ultra¬ 
violet. To check this, we have re-determined the renormalization parameters when matching 
at the larger scale pt = 5 GeV. The best match to the lattice data in this case is virtually 
undistinguishable from figure [TJ The renormalization parameters change considerably and, 
in particular, we now have Ma/}i ~ 0.173, which again yields eq. |66| in physical units. In 
this sense, the mass parameter Ma is invariant under the renormalization group. 
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FIG. 1. The ghost form factor ;/(p 2 ) (top) and the gluon propagator Dip 1 ) ( bottom ) in Landau gauge at 
zero temperature. Lattice data is taken from Ref. Ifl3l . and the gluon mass renoramalization parameter 
is M a « 0.48 GeV. 

It must be emphasized, however, that the determination of the absolute scales is not very 
accurate and different values may produce similar results as in figure [l] Ultimately, the scales 
should be fixed by the comparision to the string tension computed from the effective action 
for the Polyakov loop, which is currently under investigation. 
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B. Numerical procedure at finite temperature 


At finite temperature, our numerical procedure is based on eq. ( |64| . The five constants 
eq. ( f65] > entering this system are, in principle, fixed through the zero-temperature counter 
terms i]o(p c ), Z and M/\. The first two of these three parameters only adjust the overall nor¬ 
malization of the zero-temperature propagators, which in turn fix the normalization of the 
finite-temperature propagators as well. In practice, the field normalization is rarely expressed 
through the zero-temperature counter terms, but rather by the alternative condition that the 
propagators at all temperatures should pass through the same point at a very large scale, 

co ± (0,poo) = co\\{0,pco) = Zp 2 ^, ^00 >1. (67) 

This is based on the idea that finite temperature effects should become immaterial in the 
deep perturbative region where 0(4) invariance is restored. The lattice study in Ref. Bll . to 
which we compare, employs a similar prescription with poo = 5 GeV. If we adopt eq. ( |66| , 
this means p 0 0 ~ 5p, but much larger values are also accessible in our case. In any case, 
the prescription eq. ( j67| can be imposed simply by replacing p —> goo in the second and 
third eq. ( |64| . The values at the reference scale, Cdj_ (0, u ) and (0, //) are then temperature- 
dependent predicitons of the calculation, which should agree with the results of the more 
complicated approach |A4| if the scale goo is chosen large enough. This reasoning allows us 
to compare to lattice data conveniently, and to determine the two constants co±(0,p) and 
d>_ (0, p) in eq. ( |64) > self-consistently. 

For the ghost form factor, we would like to use the same procedure but this is hampered by 
the fact that we must impose a boundary condition in the deep infra-red in order to discern 
scaling and decoupling types of solution. This means that we have to guess the (temperature- 
dependent) intercept g(0,g c ) at Pc ~ 0 such that the solution passes through the common 
value 

?/(0, f/ 00 ) = Z c (68) 

for all temperatures. This leaves us with the remaining mass terms 

ca± (0, p 0 ) = Z m]_ , cDy (0, p 0 ) = Z m|, AM| = Z Am|. (69) 

which cannot be determined self-consistently by any scaling procedure, as they are not re¬ 
lated to field renormalizations. We will determine them by fitting to the lattice data, i.e. we 
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have three parameters to cover the entire momentum range for all three propagators at each 
temperature. The need to fit eq. (69) (when we actually have equations to compute them) 
is only due to numerical issues. We can, however, take the alternative point of view that 
the coefficients eq. (69) are part of our variational ansatz (which includes the renormalization 
procedure), and the fit corresponds to a solution of the corresponding gap-equation. This is 
very much in line with the general philosophy of our approach, although it would, of course, 
be desirable to compute eq. (69) self-consistently without input from the lattice. 

It should also be mentioned that the mass parameter AMjj could actually be computed from 
eq. (A5). We have still included it in the set of fit parameters, because eq. (A5) as well as 


the full dynamical system eq. (64' contains the coupling constant Ng 2 . From eq. (64' and the 
explicit form of the loop integrals, it is apparent that the factors of Ng 2 could be eliminated 
by the same scaling procedure that we used in the case of zero temperature. However, Ng 2 is 
temperature-independent, and so, only temperature-independent norm changes (rescalings) 
of the propagators are allowed. In addition, such rescalings do not leave the mass coefficients 
in eq. (69) invariant. This means that all three mass coefficients must be fitted to lattice data 
(or ideally determined from appendix |A 4) if we rescale the propagators or change Ng 2 . In 
particular, the mass coefficients eq. (69) must be re-determined if we change the colour group 
SU(N). The bottom line is that we are entitled to (i) set Ng 2 = 1, (ii) re-determine the coeffi¬ 
cients in eq. (69) for each temperature and (iii) allow for temperature-independent rescalings 
of the propagators. Since the mass parameter Amu has negligable influence on the final result, 
we will always set Am || = 0. Then, the entire system has only two free parameters, in and 
mu, at each temperature, in addition to a temperature-independent overall norm scale for each 
form factor. 


Numerically, the main issue in eq. (64) as compared to the zero-temperature case is the ap¬ 
pearance of the Matsubara sum which replaces the frequency integral. There are various 
techniques to compute this sum. A direct evaluation is hampered by the fact that only a rel¬ 
atively small number of frequencies can be included in the coupled integral equation system 
to keep the overall computational effort under control. Typically, we include up to n max ~ 20 
frequencies which corresponds, at a typical temperature }if> « 1, to a frequency cutoff of 
2nn max / f ~ 130 ja which is much smaller than the spatial momentum cutoff A /ji = O(10 5 ). 
Whether or not this is sufficient depends on the temperature itself: At high temperatures, the 
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Matsubara sum converges quickly and summing n max frequencies is sufficient. At lower tem¬ 
peratures, more frequencies contribute and we have to restort to more elaborate summation 
techniques as explained below. 

C. Cutoff independence 

Before presenting our numerical results, we must ensure that the finite-temperature system 
eq. ( f64] > is indeed cutoff-independent. We must show this for both the momentum and fre¬ 
quency cutoff independently, and for all components of the solution at various temperatures. 
This gives a bewildering number of plots and we only present two representative cases, the 
remaining ones all displaying the same behaviour. 

Since the cutoff-dependence can be studied at each temperature independently, we do not 
need to guess the correct ghost intercept // (0, }i c ) to ensure eq. ( |68) at all temperatures. Instead, 
we take an arbitrary fixed value t](Q,ji c ) = 14 at fi c = 0. The gluon normalization is fixed at 
f/oo =4 // with Z = 0.8 and for the mass parameters, we take rn_ = m = Ma ~ 0.5 // at 
}Iq = ///10 (and Am = 0). With this setup, we study the ghost form factor ?/(/co, |k|) as a 
function of the spatial momentum k = |k|, at vanishing Matsubara frequency ko = Op] We 
vary both the frequency and momentum cutoff, at various temperatures, and check how this 
affects the results. 

In table [Ij we show the dependence of the ghost form factor on the spatial cutoff A for two 
typical temperatures /3 = 0.5/// and /l = 20.0///. We have fixed the number of Matsubara 
frequencies to n max = 20 and perform the frequency sum naively without further extrapola¬ 
tion. (We will see below that this is sufficient for these temperatures.) As can be clearly seen, 
the dependence on the spatial momentum cutoff A is almost negligable. The sensitivity is 
slightly more pronounced for the lower of the two temperatures, and for larger external mo¬ 
menta k = |k|. In both cases, the maximal change of the ghost form factor rj(0,k) observed at 
k = 10// is less than 4% when varying A by three orders of magnitude. This clearly indicates 
that both the quadratic and logarithmic divergences have been removed in our formulation 
eq. (|64j), and that the perturbative region sets in fairly quickly, even below A = 10//. 

For the Matsubara sum, the situation is slightly more complicated. In general, the loop con¬ 
tributions from the Matsubara frequency #//, after integrating over angles and spatial loop 
11 Higher Matsubara frequencies ko > 0 show the same behaviour. 
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k = 0.01 

k = 0.05 

k = 0.1 

k — 0.5 

k = 1.0 

Jfc = 5.0 

k = 10.0 

A = 1 

11.03 

6.30 

4.41 

1.88 

1.51 



A = 10 

11.03 

6.30 

4.41 

1.88 

1.50 

1.25 

1.22 

A = 100 

11.03 

6.30 

4.41 

1.87 

1.50 

1.24 

1.21 

A = 10000 

11.01 

6.28 

4.41 

1.88 

1.51 

2.25 

1.22 

A = 1 

13.82 

13.09 

12.33 

9.09 

7.69 



A = 10 

13.82 

13.09 

12.33 

9.12 

7.65 

5.53 

5.00 

A = 100 

13.82 

13.09 

12.34 

9.21 

7.82 

5.72 

5.19 

A = 10000 

13.82 

13.09 

12.34 

9.20 

7.78 

5.68 

5.20 


TABLE I. Dependence of the ghost form factor t](0, k) on the spatial momentum cutoff A. The upper 
table is for [6 = 0.5, the lower one for /3 = 20.0. All dimensionfull quantities are measured in units of 
the renormalization scale }i. 




FIG. 2. The contribution A u (k) of the n'th Matsubara mode to the gluon form factor i] (0, k) at external 
momentum k. The left panel shows data for k — 0.5 }i and various temperatures on a logarithmic plot, 
while the right panel shows the same for k — 10 //. The contributions are multiplied by n 3 to better 
show their asymptotics. 

momentum, decay as 1 /n 3 asymptotically, which leads to a convergent Matsubara sum j^] To 
demonstrate this behaviour, let us consider the loop contribution A n (k) to the ghost form 


12 Our counting is such that n > 0 and contributions with n / 0 are actually the sum of n and (— n), which are 
equal in the present case. 
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factor from the first line in eq. (641, 


7 ( 0 rk)- 1 


7(°/?0 1 = 


If] ( 0 / }i-c) 


Lt]{0, He) L v (0,k) 


00 P OO /* 1 00 

= r r L dc l dz ■■■='£ A«(*), (70) 

n=0 Jo J ~ 1 n =0 

where the dots indicate the integrand of the corresponding loop integrals (cf. appendix |A 3[ >, 
and we have put the external Matsubara frequency fco = 0 for convenience. In figure [2j we 
plot A„ multiplied by n 3 , as a function of the Matsubara summation index n. As can be clearly 
seen, the asymptotic behaviour A n (k) ~ n~ 3 is reached very quickly at small momenta (left 
panel, k = 0.5 ft), for all temperatures down to /3 = 20/ ft. The same limit requires much 
more frequencies (up to 40) at large external momenta as demonstrated in the right panel 
(Jfc = 10 ft). The asymptotics of the individual contribution allows for an accurate estimate of 
the remainder in the infinite Matsubara series. However, the n = 0 term still dominates in all 
cases J^] and the partial sum is found to always saturate once the asymptotic behaviour has set 
in. For instance, the worst convergence in figure [2] is at k = 10 ft for the lowest temperature 
ft = 20/y (right panel). In this case, the asymptotics is not completely reached at n = 20, but 
the partial sum of the first 20 frequencies still agrees with the value extrapolated from the full 
asymptotic behaviour to better than 0.1%. 

In figure [3j we show again A n (k) n 3 , but this time on a linear plot which exhibits the slower 
convergence to the asymptotics more clearly. In the left panel of figure [3j we take a large mo¬ 
mentum k = 10 }i and study two different temperatures. It is apparent that n = 30 frequencies 
are enough to see the asymptotic behaviour and saturate the Matsubara sum at the intermedi¬ 
ate temperature /3 = 20/ fi. As we further lower the temperature to j 6 = 100/ ft, no sign of the 
correct asymptotics can be seen within the first 40 frequencies, although the n = 0 term still 
dominates in this case. The right panel of figure [3] shows that the convergence speed depends 
very much on the external momentum and, quite generally, decreases as the momentum is 
increased. 


The bottom line is that summing up to n = 30 frequencies directly is sufficient to saturate the 
Matsubara sum for all relevant momenta down to temperatures of about /3 = 20/ }t, though as 
few as 4 frequencies are necessary at higher temperatures. By contrast, very low temperatures 

13 This is not visible in the plot due to the factor n 3 . 
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FIG. 3. The contribution A n (k) of the n'th Matsubara mode to the gluon form factor 7/(0, k) at external 
momentum k. The left panel shows data for a large momentum k = 10 }i and two different temperatures 
on a linear plot, while the right panel displays data for a very low temperature j6 = 100 / }i at two 
different momenta k. Again, the contributions are multiplied by ft 3 to better display their asymptotics. 

can only be trusted for small momenta k <C pi, and would require an excessive amount of CPU 
time otherwise^] We will mainly restrict our numerical effort to include up to n = 20,... ,40 
frequencies at maximum, which means that we must content ourselves with temperatures not 
much below fi = 20/pi. 


D. Comparision with lattice results 

Finally, figure |4] and [5] show our cummulative results for the ghost and gluon propagator, 
respectively, at various temperatures compared to the SU( 3) lattice data taken from Ref. Ifl3l . 
The lattice data was renormalized at a large scale pica = 5 GeV in the perturbative region, with 
rj{0,}loo) = ?4= D -L(0,j«oo) = pilo D|| (0, pica) = 1. 

To translate this setting in our notation, we first set the arbitrary renormalization scale pi = 
pica = 5 GeV, because pi is anyhow redundant once we change our renormalization prescription 
to match the lattice convention. In these units, the critical temperature is at fi* = 16.67. The lat¬ 
tice data goes down to temperatures of about T = 0.67 T* , which translates into fi = 25 which, 
according to the studies in the previous section, is still treatable with reasonable numerical ef¬ 
fort. The gluon propagator normalization becomes Z = cu(0,pico)/pi 2 = 1/ (D(0,pica) pi^a) = 1- 
Finally, the intercept 7/(0, pi c ) is chosen at each temperature such that;/(0, pica) = 1 as discussed 

14 Series accelerators also do not help in this case, since the eventual asymptotics has not been reached, cf. fig.pl 
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FIG. 4. The renormalized ghost form factor at various temperatures, as a function of the spatial 
momentum k = |k| with frequency /cq = 0. Lattice data is taken from Ref. fll3H . 

in the previous section. This fixes all the normalization, and leaves us with just the two fit 
parameters and my (j 6)/}i at each temperature p*] 

As can be seen from the plots, our approach is able to reproduce the basic properties of the 
propagators in full Yang-Mills theory at finite temperatures: 

1. the gluon propator is suppressed as the temperature increases 

2. the ghost form factor is enhanced as the temperature increases 

3. the temperature sensitivity of the longitudinal gluon propagator Dm (0, k ) is much greater 
than for the transversal propagator (note the logarithmic scale in figure [5]). 

Quantitatively, the gluon propagator seems to miss some strength in the intermediate to low 
momentum regime, while the ghost form factor is too steep in this region. These problems 
seem, however, to be related to the determination of the overall scale }i = jico, which was 
taken from the match of the gluon propagator with the lattice data at pi = 5 GeV. As can be 
seen from figure [5j this determination is not very accurate as the propagators lie on top of 
each other for all temperatures at least down to about 3 GeV. Thus, there is some uncertainty 
of almost a factor of 2 in fi, which translates in the other dimensionfull quantities such as the 

15 We set Amu = 0 throughout as it has no effect on the final results. 
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norm of the gluon propagator or the mass parameters. Eventually, this scale should better be 
fixed intrinsically by computing a dimensionfull quantity such as the string tension from the 
Polyakov loop potential. 


VI. SUMMARY AND CONCLUSIONS 

In this paper, we have extended the covariant variational principle for Yang-Mills theory in 
Landau gauge to the case of non-zero temperatures. In the course of the derivation, we 
have also clarified and simplified the renormalization procedure at zero temperature, and re¬ 
determined the gluon mass parameter M,\ (which is the only free parameter in our approach 
besides the overall normalization of the propagators) by comparision to lattice data. At finite 
temperatures, no further counter terms are necessary to put the system in a form that is 
manifestly independent of the UV cutoff, both in spatial and temporal momentum direction. 
We have carefully corroborated this cutoff-independence through numerical studies. Two 
mass parameters, which in principle are fixed from the T = 0 sector, had to be fitted to lattice 
data for practical reasons. The results are in qualitative agreement with recent high precision 
lattice data. We have also discussed possible causes of the discrepancy, which cannot be fixed 
reliably without an intrinsic determination of the scale }i. 

The present framework allows for a number of further invstigations. For instance, the ques¬ 
tion of possible deconfinement phase transition in our approach should be answered from 
the effective action for the Polyakov loop, which gives a much stronger argument than the 
0(4)-symmetry violation in the propagators alone. This investigation is currently underway. 
In addition, the inclusion of fermions and a chemical potential is possible without further 
conceptional obstacles. To obtain realistic results, it may, however, be necessary to go beyond 
the Gaussian ansatz used in the present paper. This can be achieved e.g. with the Dyson- 
Schwinger method used in the context of variational calculations within the Hamiltonian 
approach in Coulomb gauge |[43l . 
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FIG. 5. The renormalized gluon propagator at various temperatures, as a function of the spa¬ 
tial momentum k — |k| with frequency Icq — 0. The left panel shows the transversal component 
D^(0,k) — d>j_(0, k)^ 1 , while the right panel shows the longitudinal component Dy (0 ,k) — cDn (0,A:) -1 . 
Lattice data is taken from Ref. m. 

Appendix A: Explicit calculations 

1. The zero-temperature gap equation 


In ref. 0 the (unrenormalized) gap equation at zero temperature was derived in the form 

cu 0 (k) = k 2 + Xo (k) + Ml Ml = Ng 2 iff 

Xo (k) = Ng 2 lff(k) ijoiky 1 = 1 -Ng 2 lff(k), (Al) 

where the index '0' on all quantities indicates that these are the zero-temperature profiles. The 
integrals in this equation are given by 

d^q 1 


A o) = 9 

M ~ 4, 


4 0) « = 
if w = 


(2tt) 4 q 2 + Ml + xo{q) 


d 4 q 

(2^F 

d A q 

w 


l-(k- q) 2 
l-(k- q) 2 


vo( k -q) 

J ( k-q) 2 u> 0 {q ) 

yo(fc-q)yo(q) 
( k~q) 2 


(A2) 


2. Zero temperature mass functions 


Here, we want to show the zero temperature limit 

lim Mj_(/3) = Mq , 

/3->oo 


lim AMj(^) = 0 

>-oo 11 


(A3) 
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of the mass functions that appear in the gap equation (441 and (471 in the main text. We 
employ the auxiliary relation 


lim I dqf(q 2 ) \ = j 

p->°o Jp uy,J q 2 4 


d^q 
(2 n) 1 - 


ftf)’ 


(A4) 


which follows from 0(4)-invariance in the limit fi —> oo. If we use this relation and take the 
zero-temperature limit eq. (52} into account, we find 

2 

Ng 2 f d A q A + B(q)_Ng 2 f d^q T — 3 


lim Mi (/!) = 
p —>oo Z 


= Ng 2 


( 2 /r) 4 cv 0 (q) 

d 4 q 1 


(2 7i) 4 cv 0 {q) 


2 + 


1 q 2 




^M 2 


(27t) 4 u>o{q) [ 3 q 2 _ 

[ A 1 

I (27t) 4 too(q) ["'34 


„ 1 3 

2 + - • - 


The second equation (A3} follows in exactly the same manner. 


3. Loop integrals at finite temperature 


The integrals in the main result eq. (62} take the explicit form 

h (k-q ) 2 a; !i( t ?) 

V(k-q ) 


L, 7 (k 0/ |k|) = dq 


p (k - q) 2 


(l-(k-q) 2 ) 0 ± (q) 1 — d>ii (q) 


-l 


li(k 0/ |k|) = - dq 


V(k-q)y(q) q 

(*-<?) 2 


2 r 


1- (k-q) 


4(*o,|k|) = dq 


q(k-q)q(q) 

(k-q) 2 


q 2 0 + q 2 (k-q)' 


- ( k ' d? 


M 2 L (fS) = 1 -Ng 2 jdq 
A M\(fS) = l -Ng 2 jdq 


4 

(2q 2 0 + 3q 2 ' 

) 1 


ti ± (q) 

v ^?o + q 2 

> ti\\(d)_ 

2 , 

^o~H 2 \ 

1 


G ± (q) 

y do + q 2 / 

G\\(q) 



(A 5) 


For simplicity, we have written q(q) etc. for the various form factors inside the loop integrals, 
although they depend, of course, on the two invariants qo and |q| separately. 
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4. Renormalization constants at finite temperature 


As stated in the main text, the renormalized integral equation system eq. ( [64] ) contains five 
unknown, temperature-dependent constants which are all determined by zero-temperature 
counter terms. In detail, these relations are 


c*>n(0,/*o) =ZMi + Mi( J 6)-M2 + M}( J 6) + 4 0) ( ?<0 )-fl(0,fio) =ZM 2 A + 


x 

r(°) / 


= /'o(jMc) * + 

= Z }l 2 -\ - 


J/(0 ,Hc) 1 = Vo(Fc) 1 - l v (0,}i c ) ~ C(Fc) ~ L v (0,Hc) 

cv ± (0, } i) = Z } i 2 + Mi (/l) - M 0 2 + 4 0) 00 - # (0/ ?0 
m||(0,fO = Z// 2 + Mi(£) -M§ + Mj(£) +4 0) (/0 -4(°'?0 = Z// 2 + --- 

mi (0, p 0 ) = Z M 2 + M 2 (/3) - M 2 + ff (;t 0 ) - (0, }i Q ) 


= zm 2 a 


(A6) 


where the dots in the last form are the finite-temperature corrections to the zero-temperature 
renormlization data. These corrections are very hard to treat numerically because the cor¬ 
responding difference of loop integrals involves the finite-temperature and zero-temperature 
form factors, respectively. 
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